﻿using System;
using System.Collections;
using System.Collections.Generic;
using System.Diagnostics;
using System.Diagnostics.Contracts;
using System.Linq;
using System.Linq.Expressions;
using System.Text;
using System.Xml.Serialization;

namespace MathService
{
    /// <summary>
    /// Полином степени N-1
    ///  a[0]+a[1]*x+a[2]*x^2+...+a[N-1]*x^(N-1)+a[N-1]*x^(N-1)
    /// где N - число элементов массива коэффициентов
    /// Нулевой элемент массива при нулевой степени члена полинома 
    /// </summary>
    [Serializable]
    //[DebuggerDisplay("GetPower = {GetPower}")]
    public class Polynom : ICloneable<Polynom>, IEquatable<Polynom>, IEnumerable<double>
    {
        /* -------------------------------------------------------------------------------------------- */



        /* -------------------------------------------------------------------------------------------- */

        public static Polynom Koefficients(params double[] a) { return new Polynom(a); }

        /// <summary>Преобразовать массив корней полинома в коэффициенты прои степенях</summary>
        /// <param name="Root">Корни полинома</param>
        /// <returns>Коэффициенты при степенях</returns>
        public static double[] GetPowers(params double[] Root)
        {
            Contract.Requires(Root != null);
            Contract.Ensures(Contract.Result<double[]>() != null);
            Contract.Ensures(Contract.Result<double[]>().Length == Root.Length + 1);

            //Массив коэффициентов при степенях
            var lv_Power = new double[Root.Length + 1];

            var lv_N = Root.Length;
            lv_Power[lv_N - 0] = 1;
            for(var n = 0; n < lv_N; n++)
            {
                var lv_POld = lv_Power[lv_N - 0];
                for(var i = 1; i < n + 2 && i < lv_N + 1; i++)
                {
                    var lv_PNew = lv_Power[lv_N - i] - lv_POld * Root[n];
                    lv_POld = lv_Power[lv_N - i];
                    lv_Power[lv_N - i] = lv_PNew;
                }
            }
            return lv_Power;
        }

        public static Polynom Roots(IEnumerable<double> roots) { return Roots(roots.ToArray()); }


        /// <summary>Получить полином из корней полинома</summary>
        /// <param name="Root">Корни полинома</param>
        /// <returns>Полином с указанными корнями</returns>
        public static Polynom Roots(params double[] Root)
        {
            Contract.Requires(Root != null);
            Contract.Ensures(Contract.Result<Polynom>() != null);

            return new Polynom(GetPowers(Root));
        }

        public static Polynom Random(int Power = 3, double Ma = 0, double Da = 1)
        {
            var rnd = new Random();
            var a = new double[Power];
            for(var i = 0; i < Power; i++)
                a[i] = (rnd.NextDouble() - .5) * Da + Ma;
            return new Polynom(a);
        }

        public static Polynom RandomPower(int MaxPower = 5, double Ma = 0, double Da = 1)
        {
            return Random(new Random().Next(MaxPower), Ma, Da);
        }

        public static Polynom RandomWithIntCoefficients(int Power, int Ma = 0, int Da = 10)
        {
            var rnd = new Random();
            var a = new double[Power];
            for(var i = 0; i < Power; i++)
                a[i] = Math.Round(rnd.Next(Da) - .5 * Da + Ma);
            return new Polynom(a);
        }

        public static Polynom RandomWithIntCoefficients_P(int MaxPower = 5, int Ma = 0, int Da = 10)
        {
            return RandomWithIntCoefficients(new Random().Next(MaxPower), Ma, Da);
        }

        public static Polynom PolynomCoefficients(params double[] Coefficients)
        {
            return new Polynom(Coefficients.FilterNullValuesFromEnd());
        }

        [Pure, DebuggerStepThrough]
        public static double GetValue(double[] A, double x)
        {
            var N = A.Length;
            if(N == 0) return 0;
            var y = A[N - 1];
            for(var i = 1; i < N; i++)
                y = y * x + A[(N - 1) - i];
            return y;
        }

        public static Expression<Func<double, double>> GetExpression(double[] A)
        {
            var N = A.Length;
            Expression y = null;
            ParameterExpression px = Expression.Parameter(typeof(double), "x");
            if(N == 0) y = Expression.Constant(0d);
            else
            {
                y = Expression.Constant(A[N - 1]);
                for(var i = 1; i < N; i++)
                    y = Expression.Add(Expression.Multiply(y, px), Expression.Constant(A[(N - 1) - i]));
            }
            return Expression.Lambda<Func<double, double>>(y, px);
        }

        [Pure, DebuggerStepThrough]
        public static Complex GetValue(double[] A, Complex z)
        {
            var N = A.Length;
            if(N == 0) return 0;
            Complex y = A[N - 1];
            for(var i = 1; i < N; i++)
                y = y * z + A[(N - 1) - i];
            return y;
        }

        public static Expression<Func<Complex, Complex>> GetExpressionComplex(double[] A)
        {
            var N = A.Length;
            Expression y = null;
            ParameterExpression px = Expression.Parameter(typeof(Complex), "z");
            if(N == 0) y = Expression.Constant(0d);
            else
            {
                y = Expression.Constant(A[N - 1]);
                for(var i = 1; i < N; i++)
                    y = Expression.Add(Expression.Multiply(y, px), Expression.Constant(A[(N - 1) - i]));
            }
            return Expression.Lambda<Func<Complex, Complex>>(y, px);
        }

        public static Complex GetValue(Complex[] Z, Complex z)
        {
            var N = Z.Length;
            if(N == 0) return 0;
            var y = Z[N - 1];
            for(var i = 1; i < N; i++)
                y = y * z + Z[(N - 1) - i];
            return y;
        }

        public static Expression<Func<Complex, Complex>> GetExpression(Complex[] A)
        {
            var N = A.Length;
            Expression y = null;
            ParameterExpression px = Expression.Parameter(typeof(Complex), "z");
            if(N == 0) y = Expression.Constant(0d);
            else
            {
                y = Expression.Constant(A[N - 1]);
                for(var i = 1; i < N; i++)
                    y = Expression.Add(Expression.Multiply(y, px), Expression.Constant(A[(N - 1) - i]));
            }
            return Expression.Lambda<Func<Complex, Complex>>(y, px);
        }

        public static double GetPolynomValue(IEnumerable<double> A, double x)
        {
            return A.FilterNullValuesFromEnd().Aggregate(
                new
                    {
                        v = 0.0,
                        xx = 1.0
                    },
                (v, a) => new
                    {
                        v = v.v + a * v.xx,
                        xx = v.xx * x
                    },
                v => v.v);
        }

        public static Complex GetPolynomValue(IEnumerable<double> A, Complex x)
        {
            return A.FilterNullValuesFromEnd().Aggregate(
                new
                {
                    v = new Complex(),
                    xx = new Complex(1)
                },
                (v, a) => new
                {
                    v = v.v + a * v.xx,
                    xx = v.xx * x
                },
                v => v.v);
        }

        /* -------------------------------------------------------------------------------------------- */

        /// <summary>Коэффициенты при степенях</summary>
        /// <remarks>a[0]+a[1]*x+a[2]*x^2+...+a[N-1]*x^(N-1)+a[N-1]*x^(N-1)</remarks>
        private readonly double[] f_a;

        /* -------------------------------------------------------------------------------------------- */

        /// <summary>
        /// Коэффициенты при степенях
        ///   a[0]+a[1]*x+a[2]*x^2+...+a[N-1]*x^(N-1)+a[N-1]*x^(N-1)
        /// </summary>
        [XmlArray(ElementName = "a")]
        public double[] Coefficients
        {
            [Pure, DebuggerStepThrough]
            get
            {
                Contract.Ensures(Contract.Result<double[]>() != null);
                Contract.Ensures(Contract.Result<double[]>().Length > 0);
                return f_a;
            }
        }

        /// <summary>Степень полинома = число коэффициентов - 1</summary>
        public int Power
        {
            [Pure, DebuggerStepThrough]
            get
            {
                Contract.Ensures(Contract.Result<int>() >= 0);
                Contract.Ensures(Contract.Result<int>() == Length - 1);
                return f_a.Length - 1;
            }
        }

        /// <summary>Длина полинома - число коэффициентов</summary>
        public int Length
        {
            [Pure, DebuggerStepThrough]
            get
            {
                Contract.Ensures(Contract.Result<int>() > 0);
                Contract.Ensures(Contract.Result<int>() == f_a.Length);
                return f_a.Length;
            }
        }

        ///<summary>
        /// Коэффициент при степени <paramref name="n"/>, где <paramref name="n"/> принадлежит [0; <see cref="Power"/>]
        /// <see cref="Power"/> = <see cref="Length"/> - 1
        /// </summary>
        ///<param name="n">Степень a[0]+a[1]*x+a[2]*x^2+...<b>+a[<paramref name="n"/>]*x^<paramref name="n"/>+</b>...+a[N-1]*x^(N-1)+a[N-1]*x^(N-1)</param>
        public double this[int n]
        {
            [Pure, DebuggerStepThrough]
            get
            {
                Contract.Requires(n >= 0);
                Contract.Requires(n < Coefficients.Length);
                return f_a[n];
            }
            [DebuggerStepThrough]
            set
            {
                Contract.Requires(n >= 0);
                Contract.Requires(n < Coefficients.Length);
                f_a[n] = value;
            }
        }

        /* -------------------------------------------------------------------------------------------- */

        /// <summary>Полином степени N, нулевой элемент массива a[0] при младшей степени x^0</summary>
        /// <param name="a">a[0]+a[1]*x+a[2]*x^2+...+a[N-1]*x^(N-1)+a[N-1]*x^(N-1)</param>
        [Pure, DebuggerStepThrough]
        public Polynom(params double[] a)
        {
            Contract.Requires(a != null);
            Contract.Requires(a.Length > 0);
            Contract.Ensures(f_a != null);
            Contract.Ensures(f_a.Length > 0);
            Contract.Ensures(Power >= 0);
            f_a = (double[])a.Clone();
        }

        [Pure, DebuggerStepThrough]
        public Polynom(IEnumerable<double> a)
            : this(a.ToArray())
        {
            Contract.Requires(a != null);
            Contract.Requires(a.Any());
            Contract.Ensures(f_a != null);
            Contract.Ensures(f_a.Length > 0);
            Contract.Ensures(Power >= 0);
        }

        [Pure, DebuggerStepThrough]
        public Polynom(IEnumerable<int> a)
            : this(a.Cast<double>())
        {
            Contract.Requires(a != null);
            Contract.Requires(a.Any());
            Contract.Ensures(f_a != null);
            Contract.Ensures(f_a.Length > 0);
            Contract.Ensures(Power >= 0);
        }

        /* -------------------------------------------------------------------------------------------- */

        /// <summary>Получить значение полинома</summary>
        /// <param name="x">Переменная</param>
        /// <returns>Значение полинома в точке x</returns>
        [Pure, DebuggerStepThrough]
        public double Value(double x) { return GetValue(f_a, x); }

        [Pure, DebuggerStepThrough]
        public Complex Value(Complex z) { return GetValue(f_a, z); }

        public Func<double, double> GetFunction() { return Value; }

        public Expression<Func<double, double>> GetExpression() { return GetExpression(f_a); }
        public Expression<Func<Complex, Complex>> GetExpressionComplex() { return GetExpressionComplex(f_a); }

        public Polynom DivideTo(Polynom Divisor, out Polynom Remainder)
        {
            double[] result;
            double[] remainder;
            Devide((double[])this, (double[])Divisor, out result, out remainder);
            var n = remainder.Length;
            while(n >= 1 && remainder[n - 1].Equals(0)) n--;
            Array.Resize(ref remainder, n);
            Remainder = new Polynom(remainder);
            return new Polynom(result);
        }

        [Copyright("", url = "http://losev-al.blogspot.ru/2012/09/blog-post_14.htm")]
        private static void Devide(double[] dividend, double[] divisor, out double[] quotient, out double[] remainder)
        {
            Contract.Requires(!dividend[dividend.Length - 1].Equals(0), "Старший член многочлена делимого не может быть 0");
            Contract.Requires(!divisor[divisor.Length - 1].Equals(0), "Старший член многочлена делителя не может быть 0");
            if(dividend[dividend.Length - 1].Equals(0)) throw new ArithmeticException("Старший член многочлена делимого не может быть 0");
            if(divisor[divisor.Length - 1].Equals(0)) throw new ArithmeticException("Старший член многочлена делителя не может быть 0");

            remainder = (double[])dividend.Clone();
            quotient = new double[remainder.Length - divisor.Length + 1];
            for(var i = 0; i < quotient.Length; i++)
            {
                var coeff = remainder[remainder.Length - i - 1] / divisor[divisor.Length - 1];
                quotient[quotient.Length - i - 1] = coeff;
                for(var j = 0; j < divisor.Length; j++)
                    remainder[remainder.Length - i - j - 1] -= coeff * divisor[divisor.Length - j - 1];
            }
        }

        public string ToMathString()
        {
            var result = new StringBuilder();
            var N = f_a.Length;
            for(var n = f_a.Length - 1; n > 0; n--)
            {
                var a = f_a[n];
                if(a.Equals(0)) continue;
                result.AppendFormat("{0}{1}{2}{3}{4}",
                    result.Length > 0 && a > 0 ? "+" : "",
                    a.Equals(1) ? "" : (a.Equals(-1) ? "-" : a.ToString()),
                    n > 0 ? "x" : "",
                    n > 1 ? "^" : "",
                    n > 1 ? n.ToString() : "");
            }
            if(N > 0 && !f_a[0].Equals(0))
                result.AppendFormat("{0}{1}",
                    f_a[0] < 0 ? "" : "+",
                    f_a[0]);


            return result.Length == 0 ? "0" : result.ToString();
        }

        public Polynom GetDifferential(int Order = 1)
        {
            Contract.Requires(Order > 0);

            var result = this;
            while(Order-- > 0 && result.Length > 0)
                result = new Polynom(result.Select((a, i) => new { a, i }).Skip(1).Select(v => v.a * v.i));
            return result;
        }

        public Polynom GetIntegral(double C = 0)
        {
            return new Polynom(Coefficients.Select((a, n) => a / (n + 1)).AppendFirst(C).FilterNullValuesFromEnd());
        }

        public Polynom GetInversed() { return new Polynom(1) / this; }

        public Polynom ScalePolynom(double c)
        {
            var result = Clone();
            var cc = 1d;
            for(int i = 1, N = result.Length; i < N; i++)
                result[i] = result[i] * (cc *= c);
            return result;
        }

        public Polynom Substite(Polynom P)
        {
            var p = P.Clone();
            var result = new Polynom(f_a[0]);
            var i = 1;
            for(; i < f_a.Length - 1; i++)
            {
                result += f_a[i] * p;
                p *= P;
            }
            if(i < f_a.Length) result += f_a[i] * p;
            return result;
        }

        /* -------------------------------------------------------------------------------------------- */

        /// <summary>Явное приведение типов полинома к делегату функции преобразования double->double</summary>
        /// <param name="P">Полином</param>
        /// <returns>Делегат функции преобразования</returns>
        [Pure, DebuggerStepThrough]
        public static implicit operator Func<double, double>(Polynom P)
        {
            Contract.Ensures(Contract.Result<Func<double, double>>() != null);
            return P.Value;
        }

        /// <summary>Оператор сложения двух полиномов</summary>
        /// <param name="P">Первое слагаемое</param>
        /// <param name="Q">Второе слагаемое</param>
        /// <returns>Сумма полиномов</returns>
        [DebuggerStepThrough, Pure]
        public static Polynom operator +(Polynom P, Polynom Q)
        {
            Contract.Requires(P != null);
            Contract.Requires(Q != null);
            Contract.Ensures(Contract.Result<Polynom>() != null);
            Contract.Ensures(Contract.Result<Polynom>().Power == Math.Max(P.Power, Q.Power));

            var NP = P.Length;
            var NQ = Q.Length;
            int N;
            Polynom a, b;
            if(NP > NQ)
            {
                a = P.Clone();
                b = Q;
                N = NQ;
            }
            else
            {
                a = Q.Clone();
                b = P;
                N = NP;
            }

            for(var n = 0; n < N; n++)
                a.f_a[n] += b.f_a[n];

            return a;
        }

        /// <summary>
        /// Оператор отрицания полинома (изменяет знак всех коэффициентов на обратной). Эквивалентно домножению полинома на -1
        /// </summary>
        /// <param name="P">Отрицаемый полином</param>
        /// <returns>Полином Q = -P</returns>
        public static Polynom operator -(Polynom P)
        {
            Contract.Requires(P != null);
            Contract.Ensures(Contract.Result<Polynom>() != null);

            return new Polynom(P.f_a.Select(a => -a));
        }

        /// <summary>Оператор вычетания полинома Q из полинома P</summary>
        /// <param name="P">Уменьшаемое</param>
        /// <param name="Q">Вычитаемое</param>
        /// <returns>Разность</returns>
        public static Polynom operator -(Polynom P, Polynom Q)
        {
            Contract.Requires(P != null);
            Contract.Requires(Q != null);
            Contract.Ensures(Contract.Result<Polynom>() != null);
            Contract.Ensures(Contract.Result<Polynom>().Power == Math.Max(P.Power, Q.Power));

            var NP = P.Length;
            var NQ = Q.Length;
            int N;
            Polynom a, b;
            if(NP > NQ)
            {
                a = P.Clone();
                b = -Q;
                N = NQ;
            }
            else
            {
                a = -Q.Clone();
                b = P;
                N = NP;
            }

            for(var n = 0; n < N; n++)
                a.f_a[n] += b.f_a[n];

            return a;
        }

        public static Polynom operator *(Polynom P, Polynom Q)
        {
            Contract.Requires(P != null);
            Contract.Requires(Q != null);
            Contract.Ensures(Contract.Result<Polynom>() != null);

            var power = P.Length + Q.Length;
            var a = new double[power + 1];
            for(var i = 0; i < P.Length; i++)
                for(var j = 0; j < Q.Length; j++)
                    a[i + j] += P[i] * Q[j];

            return new Polynom(a.FilterNullValuesFromEnd());
        }

        /// <summary>Результат деления полиномов</summary>
        public struct PolynomDevisionResult
        {
            /// <summary>Частное полиномов</summary>
            public readonly Polynom Result;

            /// <summary>Остаток деления полиновов</summary>
            public readonly Polynom Remainder;

            public readonly Polynom Divisor;

            /// <summary>Инициализация результата деления полиномов</summary>
            /// <param name="Divisor"></param>
            /// <param name="Result">Частоное</param>
            /// <param name="Remainder">Остаток от деления</param>
            public PolynomDevisionResult(Polynom Result, Polynom Remainder, Polynom Divisor)
            {
                this.Result = Result;
                this.Remainder = Remainder;
                this.Divisor = Divisor;
            }

            public double Value(double x) { return Result.Value(x) + Remainder.Value(x) / Divisor.Value(x); }

            public Func<double, double> GetFunction() { return Value; }

            public override string ToString()
            {
                return string.Format("({0}) + ({1}) / ({2})",
                    Result.ToMathString(), Remainder.ToMathString(), Divisor.ToMathString());
            }

            /// <summary>Оператор неявного преобразования результата деления полиномов в полином результата</summary>
            /// <param name="Result">Результат деления полиномов</param>
            /// <returns>Частное</returns>
            public static implicit operator Polynom(PolynomDevisionResult Result) { return Result.Result; }
        }

        /// <summary>Оператор деления двух полиномов</summary>
        /// <param name="p">Полином делимого</param>
        /// <param name="q">Полином делителя</param>
        /// <returns>Результат деления полиномов, включающий частное и остаток от деления</returns>
        public static PolynomDevisionResult operator /(Polynom p, Polynom q)
        {
            Polynom Remainder;
            var result = p.DivideTo(q, out Remainder);
            return new PolynomDevisionResult(result, Remainder, q);
        }

        public static Polynom operator *(Polynom P, double q) { return new Polynom(P.Select(a => a * q)); }
        public static Polynom operator *(double p, Polynom Q) { return new Polynom(Q.Select(a => a * p)); }
        public static Polynom operator /(Polynom P, double q) { return new Polynom(P.Select(a => a / q)); }

        /// <summary>Оператор неявного преведения типа полинома в массив вещественных значений коэффициентов</summary>
        /// <param name="p">Полином</param>
        /// <returns>Массив значений коэффициентов</returns>
        [DebuggerStepThrough]
        public static explicit operator double[](Polynom p) { return p.Coefficients; }

        [DebuggerStepThrough]
        public static implicit operator Polynom(double a) { return new Polynom(a); }
        [DebuggerStepThrough]
        public static implicit operator Polynom(float a) { return new Polynom(a); }
        [DebuggerStepThrough]
        public static implicit operator Polynom(int a) { return new Polynom(a); }
        [DebuggerStepThrough]
        public static implicit operator Polynom(short a) { return new Polynom(a); }

        #region Implementation of IEquatable<Polynom>

        public bool Equals(Polynom other)
        {
            if(other == null) return false;

            var a = f_a;
            var b = other.f_a;
            var power = a.Length;

            if(power != b.Length) return false;
            for(var i = 0; i < power; i++)
                if(Math.Abs(f_a[i] - other.f_a[i]) > double.Epsilon)
                    return false;
            return true;
        }

        public Polynom Clone() { return new Polynom(f_a); }
        object ICloneable.Clone() { return Clone(); }

        [Pure, DebuggerStepThrough]
        public override bool Equals(object obj) { return Equals(obj as Polynom); }

        #endregion

        [Pure, DebuggerStepThrough]
        public override int GetHashCode()
        {
            return f_a.Select((i, a) => i.GetHashCode() ^ a.GetHashCode()).Aggregate(0x285da41, (S, s) => S ^ s);
        }

        /// <summary>Возвращает перечислитель, который осуществляет перебор элементов коллекции.</summary>
        /// <returns>
        /// Объект <see cref="T:System.Collections.IEnumerator"/>, который может использоваться для перебора элементов коллекции.
        /// </returns>
        /// <filterpriority>2</filterpriority>
        IEnumerator IEnumerable.GetEnumerator() { return ((IEnumerable<double>)this).GetEnumerator(); }

        /// <summary>Возвращает перечислитель, выполняющий перебор элементов в коллекции.</summary>
        /// <returns>
        /// Интерфейс <see cref="T:System.Collections.Generic.IEnumerator`1"/>, который может использоваться для перебора элементов коллекции.
        /// </returns>
        /// <filterpriority>1</filterpriority>
        IEnumerator<double> IEnumerable<double>.GetEnumerator() { return f_a.Cast<double>().GetEnumerator(); }

        public override string ToString()
        {
            return f_a.Aggregate(new StringBuilder(), (S, a, i) => S.AppendFormat("{0}{1}{2}",
                                        a < 0 || i == 0 ? "" : "+",
                                        a,
                                        i == 0 ? "" : string.Format("*x{0}", i == 1 ? "" : "^" + i))).ToString();
        }
    }
}
